# Load necessary packages
library(dplyr)
library(ggplot2)
library(did)

# Load necessary data

load("00_data/01_data_processed/data_epa_analysis.Rdata")

# Preparing the data for modeling

# Set treatment year for analysis
epa_only$treatment_year_reg <- ifelse(epa_only$treatment_region == 0, 0, 1995)
epa_only$treatment_year_sta <- ifelse(epa_only$treatment_state == 0, 0, 1995)

# Combine outcome variables
epa_only$outcome_any <- 
  ifelse(
    epa_only$outcome_epa_acid_included == 1 |
      epa_only$other_inspection_acid_included == 1,
    1,
    0
  )

# Convert PGM_SYS_ID to a numeric identifier
epa_only$PGM_SYS_ID_num <- as.numeric(as.factor(epa_only$PGM_SYS_ID))

# Main analysis for EPA inspections using regional treatment year
set.seed(123)
model_epa_region_any <- att_gt(
  yname = "outcome_any",
  tname = "year",
  idname = "PGM_SYS_ID_num",
  gname = "treatment_year_reg",
  data = epa_only[epa_only$year < 2002 & epa_only$year > 1989, ],
  alp = .01
)

# Aggregate treatment effects
aggte_model_epa_region_any_dynamic <- aggte(model_epa_region_any, type = "dynamic", na.rm = TRUE)
aggte_model_epa_region_any_simple <- aggte(model_epa_region_any, type = "simple", na.rm = TRUE)

# Plotting the results
plot_epa_region_any_dynamic <- ggdid(aggte_model_epa_region_any_dynamic)

plot_epa_region_any_dynamic$data %>%
  ggplot(aes(x = year, y = att,
             ymin = (att - c * att.se),
             ymax = (att + c * att.se))) +
  geom_point() +
  geom_errorbar(width = .1) +
  geom_vline(xintercept = -1, linetype = 3) +
  geom_vline(xintercept = 4, linetype = 3) +
  geom_hline(yintercept = 0, linetype = 1) +
  ylab("Estimate, 99% conf. band") +
  xlab("Years since program start") +
  theme_bw() +
  scale_x_continuous(breaks = -10:7) +
  scale_y_continuous(breaks = seq(-.30, .08, by = 0.04)) +
  coord_cartesian(ylim = c(-.25, .08))

save(aggte_model_epa_region_any_dynamic, file = "00_data/01_data_processed/aggte_model_epa_region_any_dynamic.Rdata")

ggsave("03_figures/appendix/figure_a12.pdf", height = 3, width = 8)

rm(list = ls())


